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Abstract 

Flood Frequency Analysis is usually based on the fitting of an extreme value distribution to the series 
of local streamflow. However, when the local data series is short, frequency analysis results become 
unreliable. Regional frequency analysis is a convenient way to reduce the estimation uncertainty. In this 
work, we propose a regional Bayesian model for short record length sites. This model is less restrictive 
than the Index Flood model while preserving the formalism of "homogeneous regions" . Performance of 
the proposed model is assessed on a set of gauging stations in France. The accuracy of quantile estimates 
as a function of homogeneousness level of the pooling group is also analysed. Results indicate that the 
regional Bayesian model outperforms the Index Flood model and local estimators. Furthermore, it 
seems that working with relatively large and homogeneous regions may lead to more accurate results 
than working with smaller and highly homogeneous regions. 

Key words: Regional Frequency Analysis - Bayesian Inference - Index Flood - L-moments - Markov 
Chain Monte Carlo 



1 Introduction 

Flood frequency analysis is essential in preliminary studies to define the design flood. Methods for 
estimating design flow usually consist of fitting one of the distributions given by the extreme value 
theory to a s ample of flood events. If modelling exceedance over a threshold is of interest, a theoretical 
justification ( Fisher and Tippett . 1 19281 : iBalkema and Haan , 1974 ; Pickandsl 19751 ) exists for the use of 
the Generalized Pareto distribution (CP). 



F{x) = 1 - 1 + 



a 



-1/6 



(1) 



where 1 + £ [x — /i) ja > 0, a > 0. fj,, a and £ are the location, scale and shape parameters. This 
distribution is defined for £ ^ 0, and can be derived by continuity in the case £ = 0, corresponding to the 
Exponential case: 

x — /i N 



F(x) — 1 — exp 



(2) 



A com prehensive review of the Extreme Value Theory is given by Embrechts et al. (1997) and Colesl 
(1200 ll) . 



However, frequency analysis can lead to unreliable flood quantiles when little data is available at the 
site of interest. A convenient way to improve estimates of flood statistics is to incorporate data from 
other gauged locations in the estimation procedures. This approach is widely applied in hydrology and is 
known as Regional Flood Frequency Analysis (RF FA). One of the most popular and simple approaches 
privileged by engineers is the Index Flood method ( Dalrvmple . 1960h . The standard procedure allows: a) 
the delineation of homogeneous regions, i.e. a set of sites which behave - hydrologically and/or statistically 
- in the same way; b) the derivation of a regional flood frequency distribution; c) and the estimation of 
the parameters and quantiles at the site of interest. 

Regions are collections of gauged basins with similar site characteristics related to the flood magnitude. 
The pooled stations are not necessarily in the proximity of the site of interest. Forming homogeneous 



regions can be achieved in various ways. Regions have first been g eograp hically established. More recent 
work promoted the use of geographically non-contiguo us regions (p3urnl. Il990 : GREHYS . 19961 ). Recent 
research has defined the concept of "region of influence" (jAcreman and Wiltshire^ 19891). Other t echniques 



can be used such as Artificial Neural Networks to identify groups of stations ( Hall et al. , 2002f) 



The Index Flood model assumes that flood distributions at all sites within a region are identical, up to 
a scale factor. The Index Flood approach is not exempt from critics as its a pplication requires strong 
assumptions. One major implicit assumption, noticed by iGupta et al.l (|1994T ). is that the coefficient of 
variation of peak fl ows is to be constant across the region. This fundamental p roperty seems not to be 



verified in practice (jRobinson and Sivapalanl . 119971 ) and not physically justified (jKatz et all [2002) 



The assu mptions of the Index Floo d model need often to be relaxed to suit the observations. For this 
purpose, iGabriele and Arnelll (jl99lf) proposed a hierarchical approach to RFFA. The skewness is still 
supposed to be constant on the whole region, but the coefficient of variation and the mean annual flood can 
vary slowly from one subregion to another. However, the two authors underlined the practical difficulty 
to delineate these subregions. 

In the Index Flood model, each observation from any site within the region have the same weight. 
However, it seems not optimal as, obviously, the most precious information come from the target site. 
Indeed target data - even short - are the only one which are "really" distributed as the target site. 

We suggest here to carry out a Bayesian approach that encompasses the classical Index Flood model 
and uses the whole data in a more efficient manner. In summary, the proposed Bayesian approach 
differs from the Index Flood model as it: a) uses the at-site information in a more efficient way since 
this approach distinguishes the target site data and the regional data; b) does not impose a purely 
deterministic relationship between sites within the region. 

The main goal of this article is to test the efficiency and robustness of the developed regional Bayesian 
model when dealing with short record length series. For this purpose, classical frequency analysis i.e. 
local and traditional RFFA will be compared to the suggested regional Bayesian approach. Section [2] 
presents a brief summary of the classical Index Flood model. Relevant theoretical aspects of Bayesian 
theory are introduced and applied to flood modelling in a RFFA context in Section [3] Section |4] describes 
the data set used to illustrate the method. Section [5] describes the procedure used to elicit the prior 
distribution. Section [6] outlines the weaknesses and strengths of each approach on a typical homogeneous 
region. Finally section [7] presents an analysis of the effect of homogeneity level on quantile estimation. 



2 The Index Flood model 

The Index Flood method states that flood frequency distributions within a particular region are supposed 
to be identical when divided by a scale factor - namely the Index Flood. Mathematically, this assumption 
is expressed as: 

Q(S) = C (S) Q (R) (3) 

where Q( s > is the quantile function at site S , the Index Flood at site S and the regional 

quantile function i.e. the dimensionless quantile function valid across the homogeneous region. 

Equation ^ is the core of the model and leads to strong constraints concerning at-site distribution 
parameters. Consequently, the shape parameter is the same throughout the homogeneous region, whereas 
the location and scale parameters have simple scaling behaviour - see Appendix A. 

Equation © is supposed to be satisfied if all sites are hydrologically and/or statistically similar. There- 
fore, one of the main aspects of this approach is to identify a homogeneous region which includes the 
target site. 

Similarity in basin characteristics is necessa ry but not sufficient t o ensure the homogeneity of the region 
regarding the statistics of the flood peaks. iHosking and Wallisl (Il993l . Il997h suggested a heterogeneity 



measure Hi to assess if a region is "acceptably homogeneous" (Hi < 1), "probably heterogeneous" 
(1 < Hi < 2) or "definitively heterogeneous" (Hi > 2). Note the case Hi < seems to detect correlations 
between sites within the region. 

Once the region satisfies the homogeneity test of IHosking and Wallisl (Il993l . Il997t ). the regional flood 



frequency distribution and the related at-site distribution is computed in a classical way. That is, by fitting 



the regional distribution to the weighted mean of sample ^moments. Details f or computing heterogeneit; 
statis tics, regional flood frequency and at-site distribution can be found in ( Hosking and Wallisl , ll99S 



By definition of the Index Flood model, it can be seen that any realisations of each samples have the same 
weight. Giving equal weights to all site observations is debatable since the most relevant information is 
certainly the target site one. Relevance of the target site information is obvious as this is the only one 
which is "really" distributed as the target site. Thus, in this approach the available information is not 
efficiently used. 



3 Regional Bayesian model 



The Bayesian con cepts have already been applied with success to the regional frequency ana l ysis o f 
extreme rainfalls ( Coles and Pericchi , 20031 ) and floods ( Madsen and Rosbiere . 1997 ; Northronl 2004 ). 
Regional information is not used to build a regional distribution but to specify a kind of "suspicion" 
about the target site distribution. This is easily achieved in the Bayesian framework through the so 
called prior distribution. 



The main goal of Bayesian inference is to compute the posterior distribution. The posterior distribution 
■n(Q\x) is given by the Bayes Theorem: 



tt(0|x) 



tt (0) tt (x; 9) 
Jq tt (0) tt (xj 0) d9 
oc 7r(0)7r(x;0) 



(4) 



where is the vector of parameters of the distribution to be fitted, 8 is the space parameter, tt (x; 9) is 
the likelihood function, x is the vector of observation and tt (9) is the prior distribution. 

In theory, the posterior distribution is entirely known but is often insolvable - because of the integral. 
One of the solutions is to fix a prior model which leads to an analytical - or semi- analytical - posterior 
distri bution and which allows the posterior distribution to be computed more easily ( Parent and Bernierl . 
20031 ). Nevertheless, the most convenient way is to implement Markov Chain Monte Carlo (MCMC) 
techniques to sample the posterior distribution. This approach avoids using a purely artificial prior model 
with no theoretical and/or physical justifications. 



For our application, the likelihood function corresponds to the CP distribution as peaks over a threshold 
are of interest. From Eq. (TJJ, if the prior distribution is known, posterior distribution can be computed 
- up to a constant. The next section describes how to define the prior distribution. 



3.1 Prior Distribution 



The prior model is usually a multivariate distribution which must represent beliefs about the distribution 
of the parameters i.e. fi, a and £ prior to having any information about the data. 

As the proposed model is fully parametric, the prior distribution tt(9) is a multivariate distribution entirely 
defined by its hyper parameters. In our case study, the marginal prior distributions were supposed to 
be independent lognormal for both location and scale parameters and normal for the shape parameter. 
Thus, 

tt(0)oc Jexp[(0'- 7 ) T I]- 1 (0'-7)] (5) 

where 7, E are hyper parameters, 0' = (log^,loger, £) and J is the Jacobian of the transformation from 
0' to 0, namely J = 1/fxa. 7 = (71,72,73) is the mean vector, E is the covariance matrix. As marginal 
priors are supposed to be independent, E is a 3-3 diagonal matrix with diagonal elements (^1,^2,^3- 



3.2 Estimation of the hyper parameters 
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Hyper parameters are defined through the Index Flood concept. Consider all sites of a region except the 
target site - say the j-th site. A set of pseudo target site parameters can be computed: 

(6) 
(7) 
(8) 

for all i 7^ j, where CW is the at-site Index Flood and jjL^\a^\^ are respectively the location, scale 
and shape at-site parameter estimates from rescaled sample. 

Under the hypothesis of the Index Flood model, pseudo parameters (p,^ , , fWj for i j are expected 
to be similar to the target site distribution parameters. Note that, information from the target site sample 
is not used to elicit the prior distribution. Thus, C^' in equations ^ and {7} must be estimated without 
use of the j-th site sample. 



7i 



From these pseudo parameters, hyper parameters can be computed 

V2 =— j-^]l0gCT W , d2 = Jj—-jJ2 Var 
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(9) 
(10) 

(11) 



It is important at this step to incorporate the uncertainties on the elicitation of the prior distribution. In- 
deed, it may avoid problems related to misleading information resulting from a region not so homogeneous 
and moderating a "suspicion" that may be too true. 

For this purpose, two types of uncertainties are taken into account: the one from parameter estimation, 
and the other one from target site Index Flood estimation. Thus, hyper parameters 71 and 72 are 
estimated differently than 73 as pseudo parameters for location and scale parameters depends on the 
target site Index Flood. Under the hypothesis of independence between and |U* , ai^ the variance 
terms in Eq. ((9]) and (fT0|) are computed according these two types of uncertainties: 
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(12) 
(13) 

The independence assumption between C (J> and /u* , <rl is not too restrictive as the target site Index 

I '\ (i) (i) 

Flood C (3> is estimated independently from /xj , cr* . 



Note that Var 



of Var [logC( 



i) 



log 



(i) 



are estimated thanks to Fisher Information and the Delta method. Estimation 



is a special case and depends on the method for estimating the at-site Index Flood. 
Nevertheless, it is always possible to carry out an estimation of this variance, at least through standard 
errors. 



3.3 Specificities of the proposed prior model 



The const r uction of the prior distribution with the regional information was already suggested by 
(jNorthrod . |2004 ). Nevertheless, the location parameter - or equivalently the threshold in the GP case 
- was supposed to be known. Yet, from a theoretical point of vie w, the loca t ion p arameter can not be 
known prior to having any information from the target site sample. lNorthro"p ( 2004 ) developed a similar 
approach based on the Index Flood but uncertainty associated with the scale factor prediction was not 
considered. The prior distribution was elicited directly from the distribution of the "pseudo target site" 



estimates (/iW,6rW,fW). j n perspective, "pseudo target site" estimates are viewed as constant and 
not as random variables. When dealing with sites with a long record, uncertainties on parameter distri- 
butions are low. On the contrary, this have a much more impact for the Index Flood as uncertainties 
are as much larger as the target site Index Flood is estimated without use of at-site data - even with 
long record length sites. Note that if the prior distribution is overly accurate, esti mation and cred ibility 
intervals are influenced. For these reasons and unlike the approach proposed bv iNorthrop! (J2J304J), the 
target site Index Flood in the proposed methodology is considered to be a random variable and not a 
constant. 

Thus, our prior distribution is not too falsely "tight fit". But it reflects "real" beliefs about target site 
behaviour without any use of target site sample. 

Madsen and Rosbierd (1997) and Fill and Stedinger (1998) both presented a regional empirical Bayesian 
estimator. Both models used conjugates families for prior distributions. However, even if conjugates 
families are convenient devices, they should not only be used just because computations are easier. 
In their approaches, prior distributions are elicited with quantile regression on relevant physiographic 
characteristics. 

Our approach differs differs from the two previous empirical Bayesian approaches (i.e target site sample 
is not used to elicit priors) and respects in that way absolutely the Bayesian theory. Moreover, conjugate 
priors are not considered, but priors are suited to the data. For example, the lognormal distribution 
for both location and shape parameters is justified by a physical and theoretical lower bound as: a) 
discharge data are naturally non negative; so the location parameter should also be non negative; b) the 
scale parameter is strictly positive by definition of the GP distribution. 

This prior model is quite different from the one proposed bv I Coles and Tawn (1996) who introduced a 
lognor mal prior distribution o nly for the scale pa rameter. Note t hat it is possible to work with return 
levels ( Coles and Tawnl . 199G) or return periods (|Crowderi . Il99l instead of working with distribution 
parameters. However, regional information is suited to work directly with distribution parameters. For 
other studies, such prior models could be of interest if "suspicion" is based on return levels or return 
periods. 



4 Data description 

Streamflow data were collected at 48 gauging stations in an area reaching from the 45 th to 47 th N 
and from the 3 rd to the 8 th E. The selection of the gauging sites was initially based on the 22 regions 



into which France is divided for the implementation of the Water Framework Directive (IWasson et al 



20041 ). Seven regions cover the area under study. These regions were manually delineated taking into 
account the spatial pattern of mean annual rainfall, elevation and underlying geology. All these variables 
might influence flood generatio n processes. Therefore th is division is considered as a preliminary guide for 



pooling stations. According to lHosking and Wallid(|l997l) . pre-regions were slightly altered by identifying 



discordant sites while maximising the number of site within the region and meeting the heterogeneity 
test. Finally, a set of 14 stations was selected for this study. The heterogeneity statistic for this group is 
Hi = 0.17 < 1. Consequently, the region be considered as "acceptably homogeneous". 

The dataset includes seven tributaries to the Loire River and seven gauging stations located in the French 
part of the Rhone basin (Fig. [TJ Tab. [I}. The record length of the instantaneous discharge time series 
ranges from a minimum of 22 years to a maximum of 37 years, with a mean value of 32 years. The 
drainage areas vary from 32 to 792 km 2 . Moreover, most of the gauging stations monitored first-order 
stream catchments i.e. all but two pairs of catchments are unnested. The large majorities of flood in the 
region occur during autumn and winter and are caused by heavy liquid precipitation. 

Partial duration flood series were extracted from the time series for each station. Fig. [2] illustrates time 
series for stations U4505010, U4635010 and V3015010 and their associated thresholds. Threshold levels 
were selected to extract in average around two events per year while meeting the criteria of independence 



between floods (|Lang et all 119991) 



Three stations U4505010, U4635010 and V3015010 were of particular interest because of their extended 
record length of 37 years. The time series of those three sites are displayed in Fig. [5] 



In this case study, the scale factor was set to correspond to the 1-year return flood quantile - or equivalently 
the quantile associated with probability of non exceedance 0.5. Thus, our choice for the Index Flood 
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Figure 1: Location of the gauging stations within the studied area 



Table 1: Characteristics of the stations of the homogeneous region 



Code 


Station 


Area [km 2 ) 


X (km) 


Y (km) 


Record 


K0624510 


The Bonson river at St Marcellin 


104 


744.72 


2053.90 


1971-2003 


K0663310 


The Coise river at Larajasse 


61 


770.67 


2072.11 


1971-2003 


K0704510 


The Toranche river at St Cyr 


62.3 


752.63 


2076.68 


1977-2003 


K0813020 


The Aix river at St Germain Laval 


193 


729.48 


2093.71 


1973-2002 


K0943010 


The Rhins river at Amplcpuis 


114 


754.52 


2111.10 


1973-2003 


K0974010 


The Gand river at Neaux 


85 


743.45 


2107.75 


1972-2003 


K1004510 


The Rhodon river at Perreux 


32 


738.40 


2116.64 


1973-2003 


U4505010 


The Ardieres river at Beaujeu 


54.5 


773.67 


2130.75 


1969-2003 


U4624010 


The Azergues river at Chatillon 


336 


779.07 


2099.72 


1970-2003 


U4635010 


The Brevenne river at Sain Bel 


219 


775.90 


2092.57 


1969-2003 


U4644010 


The Azergues river at Lozanne 


792 


782.56 


2098.09 


1981-2003 


V3015010 


The Yzeron river at Craponne 


48 


785.47 


2084.50 


1969-2003 


V3114010 


The Gier river at Rive de Gier 


319 


780.54 


2062.67 


1981-2003 


V3315010 


The Valencize river at Chavanay 


36 


786.54 


2048.60 


1978-2003 



U4505010 U4635010 V3015010 
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Figure 2: Times series for sites U4505010, U4635010 and V3015010 and thresholds associated 



° I 1 1 1 1 1 ° I 1 1 1 1 1 ° I 1 1 1 1 1 1 

3.0 3.5 4.0 4.5 5.0 5.5 1 2 3 4 5 6 -0.4 0.0 0.2 0.4 0.6 0.8 

Location Scale Shape 

Figure 3: Histograms of pseudo target site estimates of location, scale and shape parameters for site 
U4505010 
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Figure 4: Regression on the basin area for estimating the at-site index flow for station U4505010 



is close to the sample m edian which was the reference in IRobson and Reedl (Il999h but differs from 
Hosking and Wailisl ( 1997 ) where the sample mean was used. This particular choice for the Index Flood 



is not unintentional as estimating the quantile with probability of non exceedance 0.5 is more robust than 
estimating the sample mean. Analysing the influence of Index Flood selection is beyond the scope of this 
work. The main point is to keep the same Index Flood throughout the case study to compare approaches 
on the same basis. 



5 Elicitation of the prior distribution 

To estimate the target site Index Flood, the most popular way is to develop an empirical formula that 
relates the flow statistic to geomorphological, land-use and climatic descriptors. This relationship is 
usually established by multivariate regression procedures. In our case study, we consider a simple model 
for which only one explanatory variable is introduced in the regression analysis: the drainage area. The 
power form model is adopted: 

C = aA b (14) 

where A is the area of the catchment. Parameters a and b are computed through ordinary least square 
procedures on logarithmically transformed data. 

However, more sophisticated models could be carried out. Nevertheless, for our case study, observations 
demonstrate that Eq. I|14p is a good parametrisation for estimating the Index Flood. Fig.[3]and[l]illustrate 
the efficiency of the regressive and prior model for site U4505010 for which: 

6 = O.m 101 (i? 2 = 0.94) (15) 



Regional information was incorporated in the prior distribution through the Index Flood model. More- 
over, uncertainties in the prior distribution were incorporated. Thus, the prior information is, on one 
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Figure 5: Proposed prior (dashed line), proposed posterior (solid line) and posterior from an uninformative 
prior (dotted line) marginal densities for GP parameters. Site U45050f0 with 5 years record length. 
Vertical lines denotes benchmark values. 

hand, not too falsely accurate and on the other hand, informative enough because of the supposed ho- 
mogeneity of stations. 



6 Performance of the Bayesian model on a homogeneous region 

When making classical inference on small samples, the uncertainties may be too large. If an extremal 
event or too many "regular events" in this short record period are present, estimation can be impacted. It 
could lead to a dramatic overestimation or underestimation of quantiles corresponding to different return 
levels. A perfect model is expected : a) to perform well enough even with small samples; b) to be robust 
enough when an extreme event occurs in the sample; c) to be robust enough when too many "regular" 
events occur in the sample. 

In this section, three different models will be applied. For this purpose, the three stations U4505010, 
U4635010 and V3015010 were selected to assess the robustness and efficiency of the local, regional and 
Bayesian regional models. These three different approaches correspond to : a) local: fit the GP distribu- 
tion to the peaks over threshold data with the Maximum Likelihood Estimator (MLE), Unbiased Prob- 
ability Weighted Moments (PWU) and the Biased Probability Weighted Moments (PWB); b) regional 
(REG): fit a regional GP distribution as described in section [2] and obtain the target site distribution; c) 
regional Bayesian (BAY): elicit the prior density from regional information, then compute the posterior 
density through MCMC techniques. As an illustration of MCMC output, Fig. O displays the prior and 
posterior marginal densities for the GP parameters of the proposed model. Marginal posterior distribu- 
tion obtained from an uninformative prior model arc also displayed. That is with the same prior model 
but with a large variance - i.e. df — 1000, i = 1 ... 3. Fig. [5] shows the relevance of regional information 
as the proposed prior model is clearly more accurate than an analysis directly from data. Moreover, 
for the proposed model and even with only 5 years record length, marginal posterior densities are more 
accurate than marginal prior densities - except for the shape parameter. Thus, combination of regional 
and target site information at two different stages is worthwhile, even when only few data are available. 
Location parameter is a special case as the modes of both marginal prior and posterior densities seem to 
be significantly dissimilar. 

As the main goal of this work is to compare models on small samples, efficiency will be evaluated on 
sub-samples from the original data. Local Maximum Likelihood Estimation on the whole sample will be 
used as benchmark to assess the performance of each model. This particular case will be denoted THEO 
in the following sections. The choice of MLE estimate as a benchmark value is reasonable because of 
its theoretical motivation and asymptotic efficiency. Moreover, the MLE approach allows the calculation 
of profile confidence intervals. This is a key point as these profile con fidence interv als are often more 
accurate than those based on the Delta Method and Fisher Information (IColesLl200lh . 



Furthermore, as interpretation on quantile estimates is more natural than for distribution parameter 
estimates, the analysis will focus on quantiles corresponding to return period 2, 5, 10 and 20 years. 
Benchmark values for these quantiles - and their associated 90% profile likelihood confidence intervals 



Table 2: Benchmark values for 2, 5, 10 and 20 years quantiles and the associated 90% profile likelihood 
confidence intervals in bracket 



Station 


Q 2 


Q 5 




Q20 


U4505010 
U4635010 
V3015010 


10.8(10.1, 11.7) 
33.0(30.0, 36.5) 
7.5 ( 6.9, 8.3) 


15.3 (13.9, 17.4) 
52.2 (45.5, 62.5) 
11.7 (10.4, 13.7) 


19.5(17.2, 23.4) 
72.2(60.2, 95.4) 
15.9(13.6, 19.9) 


24.4(20.6, 31.5) 
98.9(69.2, 200.5) 
21.3(17.3, 28.8) 
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Figure 6: Comparison of frequency curves for site V3015010 with 15 last years recording 



are detailed in Tab. [2j Benchmark values with return periods greater than 20 years will be considered 
unreliable - as uncertainties on these quantiles are too large with only 37 years of record. 

Moreover, for such return periods, benchmark values are quite equivalent to those obtained with PWM 
estimates - with a mean bias of 0.89%. So, performance of each model is not too much impacted by the 
choice of the MLE estimator for benchmark values. 

Different frequency curves for site V3015010 with only the last 15 years recording are displayed in Fig. [6] 
Let us focus on the largest observation. Return period related to this event is very high for the REG 
approach. All other models lead to significantly lower return periods. This flood event is extreme at 
regional scale but not anymore in a local context. This underestimation is due to the misuse of the 
target site sample to establish the regional distribution. On the other side, the regional Bayesian model 
performs well for all return periods. Indeed, Fig. [6] indicates that the return level curve is very similar to 
benchmark one. This is quite logical as it adds up the advantage of using efficiently the target site sample 
and a good "suspicion" on the global behaviour of the flood peak distribution thanks to the so-called 
prior distribution. Local approaches suggest a very heavy tail as the extremal event of year 2004 (see 
Fig. [51 right panel) was in the last sequence of 15 years of records. 

As one of the main goals of a RFFA procedure is to deal with small samples, the target site sample 
was truncated to obtain shorten periods of records of m years, m € {5, 10, 15, 20, 25, 30, 37}. Robustness 
and efficiency of the methods to converge to the parameters of the target site distribution are measured. 
For this purpose, quantile estimates corresponding to return period 2, 5, 10, 20 years - corresponding 
to non-exceedance probabilities 0.75, 0.9, 0.95 and 0.975 respectively - are pointed. The evolution of 




Figure 7: Evolution of Q2, Q5, Q10, Q20 estimates as the size increases for the site U4635010 and 90% 
profile likelihood confidence interval for the benchmark values - light blue area 



quantile estimates as a function of the record length period is presented in Fig. [7) The figure is achieved 
considering only the first m years ; that is, for example, estimates related to the 5-year record length 
corresponds to the period 1969 - 1973. 

Because of the extreme event observed in 1983 (see Fig. [3J middle panel), systematic underestimation of 
benchmark values for local and REG approaches can be noticed. This result shows that: on one hand, 
for small samples classical inference like MLE, PWB and PWU are too responsive if too many "regular" 
events occurred. On the other hand, for the Index Flood model, underestimation of quantiles is related 
to the underestimation of the scale factor in Eq. ([3]) because of these "regular" events. Only the 
Bayesian model performs well enough even with record lengths lower than 15 years. 

A monotonic increase of the design flood estimates with the sample length can be noticed in Fig. [7J 
This behaviour is easily explained by Fig. [51 middle panel. Indeed, only the last part of the time series 
shows really extreme events. As the record length increases, much more extreme events occur leading to 
higher estimates. The Bayesian approach is the only one which does not really present this monotonic 
behaviour. 

Moreover, the Bayesian approach is by far the most robust and accurate model as, on the whole range 
of record length, and for all benchmark values, estimation lies in the 90% profile likelihood confidence 
interval. This is not true with any other model. The advantage of incorporating regional information 
within a Bayesian framework is certainly to define a "restricted space" where distribution parameters 
belong to. Thus, the impact of a very extremal event - or conversely too many low-level events - should 
be regarded as an extreme event related to this "restricted space" . 

The gain of accuracy in the target site from using regional information is clearly established in section [5] 
(Fig. [6] and [7]) . The Bayesian approach seems to be robust even with small samples while being accurate 
with larger sample. The poor performance of the REG model is related to a bad selection of sites within 
the "homogeneous" region being considered and estimates may be more accurate if "better" regions were 
considered. Unfortunately, building up such region is difficult because of the purely deterministic relation 
©. As the Bayesian approach relaxes the REG model, the search for more homogeneous regions could 



Tabic 3: Heterogeneity statistics for the four region considered - statistics in bracket are obtained with 
the scale factor taken to be the 1-year quantile corresponding to non-cxccedance probability 0.5 



Region 


He+ 


He 


Ho 


Ho+ 


Hi 


7.11 (6.83) 


1.35 (1.37) 


0.17 ( 0.08) 


-0.60 (-0.67) 


H 2 


3.46 (3.38) 


1.00 (1.03) 


0.41 ( 0.33) 


-1.28 (-1.31) 


H 3 


1.40 (1.45) 


0.30 (0.28) 


-0.09 (-0.14) 


-1.14 (-1.18) 



be ineffective. The goal of the next section is to measure the potential gain, for the Bayesian model, 
against homogeneity property. 



7 Effect of heterogeneity degree on quantile estimation 

As indicated in the previous section, we focus now on the impact of the level of homogeneousness of the 
region. For this purpose, we consider four different regions - denoted He + , He, Ho and Ho + - which 
correspond to increasingly homogeneous regions according to the test of iHosking and Wallisl(|l997h . The 



Ho region corresponds to the region analysed in the previous section and described in Tab.[TJ All regions 
have 14 site except for the most homogeneous one Ho + which contains only 8 stations. He and He + 
regions are derived form Ho. One to five sites are withdrawn and replaced by other stations to obtain 
larger heterogeneity measure. The Ho + region is a sub-region of Ho. Heterogeneity statistics for these 
regions are summarised in Tab.[3l 

To evaluate the influence of homogeneousness level of a region on quantile estimation, models are assessed 
using two performance criteria : the Normalised Bias (NBIAS) and the Normalised Root Mean Squared 
Error (NRMSE). These indices are defined as follows : 

NBIAS = IV^l— ^ (16) 



NRMSE 



where k is the number of estimates of Q and Qi is the i-th estimate of the benchmark value Q. To compute 
these two indices, we fit all models on all trimmed periods of size m years - m € {5, 10, 15, 20, 25, 30}. 
Moreover, the overall performance of each model is evaluated u sing a rank score. This technique was 



already used to compare different models in lShu and Burn! (|2004r ). 



To calculate the rank score, the p models are ordered thanks to their performance indices - 1 corresponding 
to the best model and p to the worst. For each model, the scores for the different criteria are summed to 
obtain the overall rank score R for the model. For convenience, the overall rank score R is standardised 
in such a way that in lies within the interval [0, 1]: 

R s = (18) 
pq-q 

where p is the number of models being considered, and q the number of indices. A standardised rank 
score close to 1 -resp. - is associated to a model with a good - resp. poor - performance. 

Three quantiles are of particular interest Q§,Qia and Q20 - *- e - associated to probability of non- 
exceedance 0.9, 0.95 and 0.975 respectively. NRMSE, NBIAS and the standardised rank score for 
station U4635010 and a record length of 5 years are illustrated in Tab. |H Notations for different models 
in this table consist of one lowercase letter referring to the Bayesian approach b or Regional Index Flood 
r and the denomination of the homogeneity degree of the region. Only the MLE model does not use 
these notations as it is completely independent of the homogeneity level. 

Results from Tab.2]demonstrate that the Bayesian model performs quite well independently of the region 
being considered. However, this model seems to perform even better when applied to a "acceptably 
homogeneous" or "probably heterogeneous" region. For the Ho + region, the Bayesian approach performs 
poorly. This may be explained by the fact that the prior distribution is too informative and probably 



Table 4: Estimation of N RAISE and NBIAS for station U4635010 with a record length of 5 years 



Model 




N RAISE 






NBIAS 




Rank Score 


Q 5 


QlQ 


Q20 


Q5 


Q10 


Q20 


MLE 


0.33 


0.34 


0.39 


0.01 


-0.09 


-0.18 


0.26 


ITT 4- 

bHe^ 


0.16 


0.13 


0.18 


0.09 


—0.02 


—0.13 


0.65 


rHe+ 


0.27 


0.30 


0.37 


-0.12 


-0.22 


-0.31 


0.18 


bHe 


0.10 


0.07 


0.11 


0.08 


0.00 


-0.09 


0.85 


rHe 


0.27 


0.26 


0.28 


-0.03 


-0.10 


-0.17 


0.43 


bHo 


0.14 


0.09 


0.08 


0.12 


0.05 


-0.02 


0.76 


rHo 


0.27 


0.26 


0.27 


0.01 


-0.06 


-0.12 


0.58 


bHo+ 


0.29 


0.28 


0.25 


0.29 


0.27 


0.25 


0.19 


rHo + 


0.28 


0.27 


0.26 


0.02 


-0.01 


-0.04 


0.60 




Figure 8: Score evolution as a function of record length for Station V3015010. (a) REG scores, (b) BAY 
scores 



not consistent with the target site sample. This comment is yet not discrepant with the good overall 
performance of the REG model on this region. Indeed, as the prior distribution is elicited using equations 
©-([H]), and the scale factor C^> is estimated without any use of the target site sample, this can lead 
to a misleading prior distribution while the REG model performs well. The bad estimation of the scale 
factor is less important with a more heterogeneous region as the prior information is less informative, 
thus the Bayesian model performance is not highly impacted. 

On the other side, the overall rank score of the REG model increases with the homogeneity degree of the 
region. Yet, the overall rank score for the REG model never exceeds the value of 0.6 - reached for the 
Ho + region. This value remains much lower than the best rank score for the Bayesian model - i.e. 0.85. 
These results corroborate the superiority of the Bayesian approach. 

From Tab. |4l two conclusions can be established. On one hand, for small samples, the Bayesian approach 
is the most competitive model. On the other hand, results seem to indicate that there is no need to 
keep increasing the homogeneousness of the region as it increases the risk of being too confident in the 
"homogeneous region" without increasing significantly the efficiency of the model. 

These results are in line with similar results obtained for stations U4505010 and V3015010, except for 
the bad behaviour of the Bayesian model on the Ho + region. Indeed, for the other stations, the Bayesian 
model remains more efficient than the REG model within the Ho + region. However, its overall rank 
score remains stable through out the different region - He, Ho and Ho + . The "risk" to deal with too 
much homogeneous region - as Ho + - is also corroborated as the overall rank score for the Index Flood 
model for station U4505010 decreases dramatically until 0.06. Thus, the Index Flood for the Ho + region 
performs quite well for stations U4635010 and V3015010, while very surprisingly badly with U4505010. 

In Fig. El the evolution of the overall rank score as a function of the record length is illustrated for station 
V3015010. The left panel corresponds to the REG part, while the right one stands for the Bayesian 
approach. In both panel, the MLE score is also presented. Fig. [8] indicates that the evolution of the 
overall rank score is more stable for regional models - that is REG and Bayesian models - than for the 




Scale Shape 



Figure 9: Effect of bad estimation of target site Index Flood on marginal prior and posterior densities. 
Site U4635010 with 10 years record length. 



MLE. Furthermore, the benefit of increasing the homogeneity degree of the region is more relevant for 
the REG model than for the Bayesian model. Nevertheless, the worst Bayesian rank score is always quite 
close to the best REG rank score. This seems to indicate the superiority of the Bayesian approach. This 
last point is corroborated with the results corresponding to stations U4505010 and U4635010 except for 
the bHo + model for station U4635010 because of the bad estimation of the scale factor CW - as denoted 
earlier. The effect of bad estimation of the target site Index Flood on prior and thereby on posterior 
distributions is depicted in Fig. [5] From Fig. [51 it is overwhelming that the prior model is not appropriate 
- particularly for the location parameter. Prior for the shape parameter is not too false as it does not 
depend on the target site Index Flood estimate. 

As the record length increases, the MLE model becomes more and more efficient. In particular, for 
record lengths greater than 15 years, it is more effective than rHe + , rHe and rHo models. On one hand, 
for record lengths smaller than 15 years, MLE is always less efficient than Bayesian approaches and even 
significantly for bHe, bHo and bHo + models. This is quite logical as Bayesian estimation can be looked 
at as a restrictive maximum likelihood estimator - restriction being defined by the prior distribution. So, 
under the hypothesis that the prior distribution is well-defined, the "restrictive estimator" is unbiased 
and has a smaller variance. On the other hand, for record lengths greater than 15 years, MLE, bHe and 
bHo seems to be equivalent. 



8 Conclusion 



A framework to perform a regional Bayesian frequency analysis for partially gauged stations is presented. 
The proposed model has the advantage of being less restrictive than the most widely used regional model, 
that is the Index Flood. Several case studies from French sites were analysed to illustrate the superiority 
of the Bayesian approach in comparison to the traditional Index Flood and to local approaches. The 
influence of the homogeneousness level of the pooling group on quantile estimates was also considered. 
Results demonstrate that working with quite large and homogeneous regions rather than small and 
strongly homogeneous regions is more efficient. Further work can focus on the regional estimation of 
other characteristics of the flood hydrograph. For instance, a regional Bayesian model can focus on Flood 
Duration Frequency. 



All statistical analysis was carried out in the R Development Core Team ( 20051) framework. For this 



purpose, two packages were contributed to this software under the framework of the present research 
work. These two packages integrate the tools that were developed to carry out the modelling effort 
presented in this paper. The first one POT performs statistical inference on Peaks Over Thresholds, 
while the second one, RFA, contains several tools to carry out a Regional Frequency Analysis. These 
two packages are available, free of charge, at the web site http://www.R-project.org, section CRAN, 
Packages. 
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Appendix A. Properties of the Index Flood on GP parameters 

We provide in this appendix the proof for the following theorem: 

Theorem 1. Let X be a random variable GP distributed. So X has the Cumulative Distribution Function 
defined by: 



Let Y = CX where C € R+. Then, Y is also GP distributed with parameters (C/U, Ccr, £). 

Proof. Let X be a r.v. GP distributed with parameters (/j,, a, £) and Y = CX where C G R+. Then: 



So, Y is also GP distributed with parameters (fj,C, aC, £). The proof for the GEV case can be established 
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